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Abstract 

We present 'empirical' models (pressure vs. density) of Saturn's interior constrained by the grav- 
itational coefficients J2, J4, and Jq for different assumed rotation rates of the planet. The empirical 
pressure-density profile is interpreted in terms of a hydrogen and helium physical equation of state to 
deduce the hydrogen to helium ratio in Saturn and to constrain the depth dependence of helium and 
heavy element abundances. The planet's internal structure (pressure vs. density) and composition are 
found to be insensitive to the assumed rotation rate for periods between 10h:32m:35s and 10h:41m:35s. 
We find that helium is depleted in the upper envelope, while in the high pressure region (P > 1 Mbar) 
either the helium abundance or the concentration of heavier elements is significantly enhanced. Taking 
the ratio of hydrogen to helium in Saturn to be solar, we find that the maximum mass of heavy elements 
in Saturn's interior ranges from ~ 6 to 20 M®. 

The empirical models of Saturn's interior yield a moment of inertia factor varying from 0.22271 to 0.22599 
for rotation periods between 10h:32m:35s and 10h:41m:35s, respectively. A long-term precession rate of 
about 0.754" yr^^ is found to bo consistent with the derived moment of inertia values and assumed rota- 
tion rates over the entire range of investigated rotation rates. This suggests that the long-term precession 
period of Saturn is somewhat shorter than the generally assumed value of 1.77x10® years inferred from 
modeling and observations. 
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1 Introduction 



Models of Saturn's interior based on pre-Cassini values of the planet's gravitational coefficients and the 



equation of state of Saumon et al. (1995) have been presented by several authors (e.g., Fortney & Hubbard 



2003 Saumon & Guillot 20041. Nevertheless, there is still considerable uncertainty in Saturn's internal 



structure due to incompleteness and lack of precision in our knowledge of Saturn's gravitational field, the 
absence of information on Saturn's abundance of helium and heavy elements, and the uncertainty in the 
equation of state of hydrogen-helium mixtures at high pressures and temperatures. For example, the size 
and mass of Saturn's heavy element core and the depth of the transition from molecular to metallic hydrogen 
are unknown. An additional source of uncertainty in determining Saturn's internal structure has arisen with 
the realization that we do not know the rotation rate of the deep interior ( Gurnett et al.[ 2007). 
In this paper we seek to provide reference radial profiles of pressure (p) and density (p) in Saturn's interior 



based on mass, radius, gravitational coefficients from the analysis of Cassini ( Jacobson et al. 2006 1 and other 
data. These interior models are non-unique because of Saturn's unknown rotation rate and uncertainties 
in the gravitational coefficients, but as we will show, they provide tight constraints on the real distribution 
of pressure and density inside Saturn, on properties of Saturn such as its moment of inertia and precession 
rate, and on the composition of the planet. 

Because we represent the radial profile of density in Saturn as a polynomial function of the radial coordinate 
s, the mean radius in the interior, the profile is independent of uncertainties and assumptions about the 
equation of state of hydrogen-helium mixtures. The same is true about the radial profile of pressure which 
follows from integration of the hydrostatic equation using the radial density profile. Elimination of s between 
p(s) and p{s) gives an 'empirical' equation of state (EOS), p = p{p). The empirical EOS is dependent only 
on the assumed internal rotation rate of Saturn and it is uncertain only by the truncation of the gravitational 
field representation and the errors in the gravitational coefficients. However, the representation of the radial 
profile of density by a 6'th degree polynomial function of mean radius may be inadequate to account for 
all the features of Saturn's actual interior such as a density discontinuity at the surface of a heavy element 
core. With the exception of the small effect of the unknown Saturnian rotation rate, the empirical EOS we 
derive provides a unique polynomial function model of the radial distribution of pressure and density inside 
Saturn. 

In the next section interior models (p(s), pis), p{p)) of Saturn are derived using the 'theory of figures' 
(Zharkov & Trubitsyn, 1978). The interior models fit the atmospheric model of Lodders & Fegley (1998) 
and the Saturnian gravitational moments. The models make no assumptions about the planet's composition 
or its radial dependence. Since the rotation rate of Saturn's deep interior is still unknown we follow [Anderson] 



& Schubert (2007) and assume rotation periods between 10h:32m:35s and 10h:41m:35s. In section 3, the 



physical equation of state of Saumon et al. (1995) is used to infer the hydrogen to helium ratio of the planet 
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and its dependence on radius based on comparison with the empirical EOS of this paper. Section 4 presents 
the values of moment of inertia and the precession period of Saturn's rotation axis predicted by the empirical 
EOS for the range of rotation periods studied. We conclude with a general discussion of the results. 



2 Interior Models: Finding Radial Profiles of Density and Pres- 
sure in Saturn's Interior - An Empirical Equation of State 

The calculation of the interior models proceeds in three steps. First, the measured gravitational field and 
polar radius are used to obtain the reference geoid or the effective gravitational potential function U, where 

U ^ V + Q 

CM ( ^ /nA2n 

V 




E(^) >/2„P2„(cOS^) . (1) 



In (1), 1/ is the gravitational potential, Q is the centrifugal potential, and lo is the angular velocity of rotation. 
We take the rotation to be that of a solid body with constant angular velocity. Though the atmosphere is 
diS^erentially rotating, our models remain relevant as long as the differential rotation is shallow and does 
not affect the deep interior. Further in (1), (r, 0,0) are spherical polar coordinates, G is the gravitational 
constant and M is the total planetary mass. For a rotating fluid in hydrostatic equilibrium, only the even 
zonal harmonics are stimulated, and the gravitational potential V can be represented as an expansion in even 



Legendre polynomials (Kaula 1968 Zharkov & Trubitsyn 1978). For pure rotation, the longitudinal 



angle </> does not enter in U . The constants that define V for a particular planet are the gravitational constant 
times the total planetary mass GM, an equatorial radius a, and the harmonic coefficients J^m which can 
be inferred from Doppler tracking data of a spacecraft in the planet's vicinity, such as the Cassini orbiter of 
Saturn. 

The reference geoid is defined as the surface of constant effective potential V with occultation polar 
radius of 54,438 km. In the absence of any published Cassini data on occultation radii, we used the published 
Voyager value for the polar radius of the 100 mbar isosurface, along with its standard error (Lindal et al.. 



1985; Nicholson et al. 1995). The polar radius is expected to be relatively independent of rotation rate and 
atmospheric winds, and is held fixed in our models for the reference geoids. For each rotation period there is 
a different surface shape of equal gravitational potential in which the equatorial radius changes with respect 
to the rotation period. Occultation data can define the equatorial radius of the 100 mbar isosurface only 
for an assumed rotation period of Saturn. As shown by Anderson & Schubert (2007), the altitude above 
the reference geoid of the measured 100 mbar surface can vary by hundreds of kilometers at the equator, 
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depending on the assumed rotation period (Lindal et al., 1985; Hubbard et al. 19971. The mean radius R 
for the 100 mbar isosurface is defined by the radius of a sphere that has the same mean density po as Saturn. 
The mean density is dependent on the rotation rate as shown in Table 1. For a more rapidly rotating planet, 
the equatorial bulge is greater and the mean density is lower for the same total mass. 
[Table. 1] 

From the Cassini mission, the GM for Saturn is 37931208 km'^s"^ (Jacobson et al., 2006), and with a modern 



value for the gravitational constant G of 6.674215 x 10"" m^s^^j^g-i ( jGundlach fc Merkowitzf |2000| , the 
total mass of Saturn M is 5.683246 x lO'^® kg, with a fractional uncertainty in G equal to 14 ppm. The 
mean density po is defined as this total mass divided by the volume of the fifth order reference geoid. The 
smallness parameter m follows from its definition m = uj^R^/GAI (Zharkov & Trubitsyn, 1978), where uj 
is the angular velocity associated with the periods of rotation in Table 1. The similar smallness parameter 
q, given by ui'^a^ /GM, is used for the calculation of the reference geoid. Self-consistent values of a and q 
are found by iterating the calculation of the reference geoid as discussed in the caption of Table 1. The 
gravitational coefficients J2, J4 and Jg are simultaneously determined from the iterations using the values 
of the observed gravitational coefficients J2, J4, and Jq for the reference equatorial radius of 60,330 km 
(Jacobson et al., 2006) according to the procedure discussed in the caption of Table 1. Values of Jg and Jio 
for the fifth order calculation are obtained by extrapolation from the lower-degree coefficients. 

The characteristic pressure in the interior is defined by po—GMpo/R (Zharkov & Trubitsyn, 1978). Given 
the parameters of Table 1, the normalized mean radius f3 is defined by s/ R, where s is the mean radius in the 
interior, the normalized mean density ri{(3) is p{s)/po, and the normalized pressure ^(/3) is p{s)/pQ. From the 
theory of figures (Zharkov & Trubitsyn, 1978), a particular interior model is characterized by an assumed 
density distribution ri((3) and smallness parameter to, and by the shape of level surfaces over the interval 
< /3 < 1. The basic idea of the interior calculations is to start with a best guess for r]{f3), compute the level 
surfaces in the interior, and then evaluate the harmonic coefficients J2, J4, and Jg at the surface for /3 equal 
to unity. The differences between the calculated coefficients and the observed surface values from Table 1 
are used to correct the density function, and the process is iterated to convergence. Further discussion of 
the density distribution 7y(/3) and the interior model is given below. 

Once a density distribution that matches the observed gravitational coefficients is available, the pressure 
in the interior is obtained by integration of the equations of hydrostatic equilibrium and mass continuity. 
The equations to the first order in to are (Zharkov & Trubitsyn, 1978), 

1^3,,-, ,3, 

where the parameter a is the normalized mass M{(3)/M internal to a level surface labeled by the mean 



4 



fractional radius /3. The normalized axial moment of inertia for the planet 7 is also available from the 
density distribution by the integration 



where C is the axial moment of inertia. Because the internal density distribution matches the observed 
gravitational coefficients, the normalized moment of inertia from Eq. 4 is consistent with those coefficients, 
and with the rotation period. 

Atmospheric boundeiry condition 

We represent the internal density distribution by a single sixth degree polynomial with the first degree 
term missing. Such a polynomial contains as many unknowns as the degrees of freedom imposed by the 
atmospheric density profile and its connection to the interior. The density distribution near the surface is 
based on two degrees of freedom. The data for the interior consist of the three gravitational harmonics 
J2, J4, Je and the measured mass and mean radius of Saturn. The rotation rate is a free parameter, and 
with an assumed density distribution as a function of radius, there are three observational constraints on 
the polynomial. The implementation of the method of level surfaces finds a one-to-one match between the 
polynomial and the three gravitational harmonics at a given rotation rate. The sixth degree polynomial, 
with the first degree term set to zero, has six coefficients that can be fit to the five measurements (three 
harmonics plus two atmospheric constraints) by the method of nonlinear least squares. The indeterminacy 
is eliminated by imposing a condition on the polynomial that all the measured mass is included between 
the center and the surface (0 < /? < 1). Hence there are five measurements and five free parameters in 
the fitting model. Each of our models for a given rotation rate is unique. For polynomials of degree n > 6 
there are n — 6 degrees of indeterminacy in the fitting process. A steep increase of density at the center 
could be imposed as a boundary condition on a higher-degree polynomial, but the resulting best-fit model 
would then depend on that assumed boundary condition. With the sixth degree polynomial, an interior 
model that fits all the available data and that is also free of any physical constraints on the interior can be 
obtained. The derivative of the density goes to zero at the center, and it must be negative everywhere else 
on the interval < /3 < 1. Otherwise the density would decrease with increasing depth, which is a physical 
impossibility. We further eliminate the sixth degree polynomial coefficient by insisting that all the mass be 
used up during the integration of equation (3) from the center to the surface. The resulting polynomial with 
five free coefficients is. 




(4) 



r) = 




(5) 
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During the calculation of the level surfaces, we check to make sure that the derivative of the polynomial 
is everywhere negative, but that is not imposed as an additional constraint. If the derivative were positive 
anywhere on the interval, we would take that as a proof that the single polynomial is an inappropriate 
approximation to the true density distribution in the interior. 

In previous models (Anderson & Schubert, 2007) the polynomial and its derivative were set to zero at 
the surface. For purposes of better matching the interior polynomial to the atmosphere, we use the model 
atmosphere in Table 9.2 of Lodders & Fegley (1998), and derive least-squares normal equations for the 
polynomial coefficients. We interpolate in the table and derive the normalized density at four values of [3 
(0.9985, 0.9990, 0.9995, 1.0), where the value 1.0 represents the 100 mbar level. For the mean radius and 
mean density associated with the 10h:32m:35s rotation period (Table 1), the corresponding atmospheric 
densities rj are (0.0002666, 0.0001719, 0.0000945, 0.0000408). These four points are fitted by least squares 
in combination with the normal equations associated with the observed gravitational coefficients J2, J4, Jq. 
The errors on the atmospheric densities are taken at 100% and the errors on the gravitational coefficients are 
given by their converged covariance matrix from the fits to Cassini and other data (Jacobson et al., 2006). 

The partial derivatives of the atmospheric density with respect to the polynomial coefficients are simply 
the coefficients in the polynomial of equation (5). These partial derivatives are collected into a 4 x 5 matrix 
A, where each row of the matrix corresponds to a particular observed atmospheric density. The atmospheric 
least-squares problem is linear, and the polynomial coefficients that correspond to the model atmosphere can 
be obtained at once. However, the rank of the A matrix is not four, but two. The atmosphere imposes two 
constraints on the polynomial, similar to the simpler but less satisfactory assumption that both the density 
and its derivative are zero at the surface. The polynomial coefficients from the atmosphere alone can be 



obtained by singular- value decomposition of the matrix A and the computation of its pseudo inverse ( Lawson 



& Hanson 19741, but the resulting polynomial is not useful for extrapolation to values of (3 less than 0.9985. 
The three observed gravitational coefficients must also be introduced for a meaningful determination of the 
polynomial. The atmospheric model serves only as a boundary condition on the deeper interior model of 
interest. 

We form a diagonal weighting matrix W with the inverse squares of the four rj data on the diagonal 



(100% error). The atmospheric normal equations can then be written as (Lawson & Hanson 1974), 



{A^WA)x = A^Wz, (6) 

where the superscript T indicates a transpose, a; is a 5 x 1 column matrix containing corrections to the 
assumed polynomial coefficients, and z is a 4 x 1 column matrix containing the corresponding residuals to 
the atmospheric values of 77. These atmospheric normal equations are combined with the normal equations 
for the gravitational coefficients, where the gravitational normal equations are designated by a subscript J. 
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The best estimate x of the corrections to the polynomial coefficients from the two data sets is, 

{A^WA + A^jWjAj)x = A^Wz + A^WjZj 



(7) 



The gravitational Aj matrix is a 3 x 5 matrix containing partial derivatives of the three gravitational 
coefficients with respect to the five polynomial coefficients. The weighting matrix Wj is the inverse of the 
3x3 covariance matrix from the data fits for the gravitational field (Jacobson et al., 2006). The 3x1 matrix 
zj contains the residuals for the observed J2, J4, Je corresponding to the current estimate of the polynomial. 
The process is iterated until it converges. The inverse of the converged matrix {A^WA + A'^WjAj) is the 
covariance matrix for the five polynomial coefficients. It can be mapped onto the polynomial for purposes 
of obtaining error bars on the density distribution in the interior. The fractional error is largest near the 
surface where the assumed error is 100%, but where the density is small. 
[Fig. 1] 

Figure 1 shows the uncertainty in the p{p) relation for a rotation period of 10h:32m:35s. The solid line 
presents the computed p{p) relation, and the dotted and dashed-dotted curves present the p{p) relation 
when the error (uncertainty in the EOS) is added and subtracted, respectively. The area between these two 
curves represents the empirical EOS derived from the interior model. As can be seen from the figure, the 
difference is largest in the lowest pressure region. This is due to the 100% error assumed on the atmospheric 
density. As the pressure increases the error decreases and the three curves overlap. A smaller error in the 
atmospheric model would lead to a more accurate p{p) relation in the low pressure regions. Because the 
residuals for the gravitational harmonics are much smaller than their standard errors, the polynomial in the 
deep interior is insensitive to the relative weighting of the gravitational and atmospheric data. Even in the 
outer atmospheric layers, a 100% error in the atmospheric density produces only a Log 2 deviation from the 
best-fit polynomial. We conclude that the empirical density function is robust in the high pressure region 
and more uncertain in the upper part of Saturn's envelope. The covariance matrix also can be mapped onto 
the pressure distribution by the integration of equation (2) and the random error in the empirical equation 
of state can be estimated for a fixed value of the period. The systematic error caused by uncertainties in the 
rotation period of about plus six and minus one minute about our preferred period of 10h:32m:35s (Anderson 
& Schubert, 2007) is given by the spread in parameters across Table 1 and Table 2. However, the empirical 
equation of state is not particularly sensitive to this relatively large systematic error in the period, let alone 
the random error. Its determination is robust. 

Calculation of level surfaces and the gravitational hcirnionics 

With a given value of the smallness parameter m and the assumed density distribution r]{l3), the level surfaces 
for constant internal potential can be evaluated and the surface harmonics can be computed from a series 
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approximation in m to the equation (Zharkov & Trubitsyn, 1978), 

Ma" J„ = - j p{r)r"Pn{cos 0)dT, (8) 

where the integration is carried out over the volume r, r is the radius, 6 is the polar angle or colatitude, and 
Pn is the Legendre polynomial of degree n. 

We have coded the level surface theory, which is given by Zharkov & Trubitsyn (1978) to the fifth order, but 
we truncate to terms in third order. From an assumed density distribution r/(/3) (equation 5) the gravitational 
harmonics J2, J4, and Jg are evaluated and the three residuals for the zj matrix become available. The Aj 
matrix is obtained by finite differencing the polynomial coefficients one at a time and by running the level 
surface code five times to obtain an estimate of the partial derivatives. The covariance matrix WJ^ for the 
observed gravitational harmonics in units of 10~^ is given by (Jacobson et al., 2006), 



J2 

wy' = J4 



J2 J4 Je 

^ 0.0746 0.6166 1.3618 ^ 

0.6166 7.6984 23.2222 (9) 

y 1.3618 23.2222 93.0205 y 

[Table. 2] 

By iterating with equation (7) the polynomial coefficients of Table 2 are obtained as a best fit to the gravi- 
tational harmonics and the atmospheric data. The converged residuals for the harmonic coefficients also are 
given in Table 2. The fit is not perfect because there is a trade-off between the fit to the atmosphere and 
the fit to the harmonics. However, the fact that the two data sets are satisfied well within their respective 
standard errors lends credibility to the interior density distribution. 

The density distribution is found to be relatively insensitive to the assumed rotation rate for periods between 
10h:32m:35s and 10h:41m:35s. Figure 2 presents the pressure-density relation from the interior models. We 
present the results for rotation periods of 10h:32m:35s and 10h:41m:35s, the shortest and longest rotation 
periods considered. The density functions are very much alike even with a nine minute difference in rotation 
period. To see more clearly the difference in the two functions. Figure 2 divides the pressure range into four 
regions. If the density-pressure relations were presented for the entire pressure range of the planet, the two 
curves would overlap. 
[Fig. 2] 

The values obtained for the normalized axial moment of inertia (see equation (4)) in each of the interior 
models are also presented in Table 2. The moment of inertia of the Saturn models depends weakly on the 
rotation rate for the range of periods considered. 
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3 Hydrogen to Helium ratio in Saturn 



Saturn is mostly a convective mixture of hydrogen and helium with a minor amount of heavy elements 



(Hubbard 1968 Guillot et al. 19941. In this section we use the empirical EOS (pressure-density profile) 



derived from the interior models to investigate the hydrogen-helium distribution in Saturn's interior. Using 
the EOS of Saumon et al. (1995), we look for the hydrogen-helium mixing ratio that produces a pressure- 
density profile similar to the one found by the interior models. Theoretical models of Saturn's interior 



suggest that the planet contains ^^10 - 30 M0 of heavy elements (Saumon & Guillot 20041, however, the 



exact amount of the material and its composition are unknown. For simplicity, when fitting the empirical 
pressure-density models we include only hydrogen and helium. In that case, the high-Z material is manifest 



as a higher helium mass fraction (Guillot 2007 BarafFe et al. 2008), so an increase of the helium mass 



fraction above the proto-Sun value can reveal the amount of high-Z material in the interior. 
To produce a physical EOS that can be compared to the p-p relation from our interior models we take an 
adiabatic EOS of a homogeneous hydrogen-helium mixture. An adiabatic EOS is justified since Saturn's 
interior is expected to be convective, with an opacity increase with increasing pressure and temperature 



(Hubbard 19731. Although convection might be suppressed by compositional gradients, condensation or 



rotation, Saturn interior models predict either no radiative zone or a very narrow one (Guillot et al. , 2004). 
An adiabat can be produced once the value of the entropy is known. The entropy of a hydrogen-helium 



mixture is given by (Saumon et al. 1995 1, 



S{p, T, X) = XSh{p, T) + YSHe{p, T) + 5„„,(p, T, X) 



(10) 



where X is the mass fraction (mass mixing ratio, mass of atoms over the total mass) of hydrogen, Sr is 
the entropy of hydrogen, Y = 1 — X and Shc are the mass fraction and entropy of helium, respectively, 
and Smix is the entropy of the mixture. The entropy of the mixture is determined using the constraint on 
the temperature at the 1 bar level in Saturn's atmosphere. The commonly used value for the Saturnian 



temperature at the 1 bar pressure level is 134.4 K (Lindal 19921, but it has been suggested by Guillot (1999) 



that this temperature could be as high as 145 K. In addition, the temperature at the pressure of 1 bar is 
obtained from radio-occultation measurements for an assumed helium to hydrogen ratio. Temperature is not 
directly observed but inferred based on an assumption about the composition of the planet's atmosphere. 



Because of the uncertainty in the surface temperature, we follow Saumon & Guillot (2004), and take Saturn's 
temperature at 1 bar to range from 130 K to 145 K. The temperature at 1 bar is significant because it 
determines the entropy, S'(p, T, X), and therefore the adiabat. Once the entropy is determined, we obtain a 
density function p{p, X, S) that can be compared to the empirical density function derived in the previous 
section. 

We compute the adiabats over Saturn's interior pressure range for hydrogen to helium mixing ratios 
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ranging from pure hydrogen to pure helium. We determine the composition which gives the best fit, in the 
least-squares sense, of the hydrogen-helium mixture EOS to the empirical p-p relation found by the interior 
models (section 2). We find that the minimum of the standard deviation function is broad and as a result 
we also report the range in composition over which the best fit value changes by ±10% (uncertainty factor). 
This range also accounts for the uncertainty in the empirical EOS. 
[Fig. 3] 

Figure 3 shows the difference between the empirical EOS and the computed density-pressure relations which 
give the best fit for rotation periods between 10h:32m:35s and 10h:41m:35s, and different temperatures at 
the 1 bar level. The dotted, solid, dashed-dot and dashed curves represent temperatures of 130, 135, 140 
and 145 K, respectively. The larger density difference between the empirical and physical p-p relations in the 
deep interior suggests that a homogeneous hydrogen-helium mixture is insufficient for describing Saturn's 
deep interior, and that the innermost region is likely composed of heavier materials, possibly in the form of 
a heavy element core. Table 3 summarizes the best fit composition value, and the composition range that 
includes the uncertainty factor for all the considered rotation periods and surface temperatures. 
[Table. 3] 

We find that faster rotation results in lighter composition, i.e., a larger hydrogen mass fraction. Higher 
surface temperatures lead to a larger mass mixing ratio of helium. The best fit values of the hydrogen mass 
fraction are found to range from X=0.82 to 0.65. In principle, the average ratio of hydrogen to helium in 



Saturn should be similar to the abundance of the proto-Sun, X'^0.725 (Bahcall et al. 19951. Fits which 
produce a hydrogen mass fraction larger than solar can be excluded as being unrealistic. We find that several 
combinations of rotation period and 1 bar temperature can provide hydrogen mass fractions which are close 
to the proto-Sun value. In some cases, the helium mass fraction is found to be larger than the proto-Sun 
value. This enrichment in helium implies the existence of heavier elements in Saturn's interior. In these cases 
we estimate the mass of heavy elements by subtracting the helium mass based on the proto-Sun composition 
from the "enriched helium" abundance; since the mass fraction of hydrogen is smaller than solar for these 
models our results give only an upper bound for the mass of heavy elements. The upper bounds on the mass 
of heavy elements are given in Table 3. The maximum mass of heavy elements accounts for the lower bound 
of the hydrogen mass fraction over the entire composition range given in column 4. The total maximum 
mass of heavy elements is found to range from ~ 6 to 20 . 



In the following section we repeat the procedure described here for different pressure regions in the 
interior. This enables a determination of whether Saturn's composition varies with depth. We find that 
the helium to hydrogen ratio in Saturn's upper atmosphere can be significantly smaller than the proto-solar 
value possibly due to sedimentation of helium towards the center. 
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3.1 The distribution of helium 

The mass mixing ratio of helium in Saturn's atmosphere is uncertain. Radio occultation measurements and 
analysis of spectra from Voyager IRIS found that helium is depleted in Saturn's atmosphere, with Y= 0.06 



± 0.05 (Conrath et al. 19841. More recent calculations using Voyager data have led to higher values ranging 



from Y = 0.18 to 0.25 (Conrath & Gautier 20001. These values are still lower than the protosolar value 
C^proto ~ 0.275). It has been suggested by several authors (Stevenson 1975 Stevenson & Salpeter 1977a|b 



Fortney & Hubbard 2003 20041 that the depletion of helium in Saturn's atmosphere is a consequence of 



helium separation from hydrogen. If helium becomes insoluble in hydrogen, it can coagulate to form helium 
droplets that settle towards the planet's center (due to larger density). Helium separation provides an 
explanation for the low helium abundance in Saturn's atmosphere, and it also offers an additional energy 



source that seems necessary to explain the long-term evolution of Saturn ( Fortney & Hubbard 2003 ) . It is 
therefore possible that the average ratio of hydrogen to helium in Saturn is similar to the proto-Sun value, 
but that the distribution of helium is not uniform throughout the interior. 

In this section we investigate whether the empirical EOS suggests a dependence of Saturn's composition 
on depth in its interior. We focus on three different pressure regions: the first represents the uppermost 
atmosphere and is defined to be from logP{Mbar) — —6 to logP{Mhar) — —4. In this pressure region 
hydrogen is in the molecular form. The second pressure region covers most of the planetary mass (99%) 
ranging from logP{Mbar) = —3 to logP{Mhar) — 1.12. This pressure region excludes the very low pressure 
region in which our interior empirical model is most uncertain. The last region is the innermost part of the 
planet, ranging from logP{Mhar) to the center of the planet, which is found to be at logP{Mhar) ^ 1.12 
(the exact value depends on the assumed rotation rate). The transition from molecular hydrogen to metallic 
hydrogen occurs at pressure of ~ 1 Mbar, and it has been suggested that in the metallic region helium is 



most insoluble (Hubbard & Dewitt 1985 Stevenson 19821. 

Again, we consider surface temperatures that range from 130 to 145 K at 1 bar, and rotation periods between 
10h:32m:35s and 10h:41m:35s. We find that the upper envelope is depleted in helium and that the helium 
mass fraction increases significantly with depth. The best fit helium mass fraction in the low pressure region 
is found to range from to 13% in agreement with observations (Conrath et al., 1984). In the pressure region 
that covers most of the planetary mass the helium mass fraction values are found to be larger than the ones 
found in the previous section when the entire pressure range of Saturn was considered. This suggests that 
the composition of Saturn is inhomogeneous, with the mass of helium or other heavy elements increasing 
with depth. In the high pressure region, more than 70% of the material is helium. However, since only 
hydrogen and helium are considered, this large helium concentration probably represents an enrichment of 
heavier elements. The presence of a solid core would reduce the helium mass fraction to lower values than 
found here. However, our results still suggest helium depletion in the upper atmosphere regardless of the 
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composition of the deep interior (see section 5 for further discussion) . Table 4 presents the compositions and 
their uncertainty factors for all the considered cases. 

While the bulk composition (see Table 3) varies significantly with the assumed surface temperature, we find 
that the composition is insensitive to the assumed surface temperature in the high pressure regions. 
[Table. 4] 



4 Precession of Saturn's Pole 



Saturn's axis is tilted to its orbital plane. The solar torque exerted on Saturn's oblate figure and on its 
equatorial satellites results in a precession of the planet's axis of rotation. Since the orbit plane of Saturn is 
not fixed in space, the precession of Saturn's pole is not constant, but changes periodically, slowly over time 
( Ward fc Hamilton] 20041. The precession rate of a planet depends on both its moment of inertia and its 
rotation rate (and the torques driving the precession of the rotation axis). For Saturn, both quantities are 
a priori unknown. Saturn's precession rate can be determined in different ways and has been computed by 



several authors. Bosh (1994) obtained the precession rate by a combination of Voyager occultation data and 
ground-based stellar occultation (28 Sgr), suggesting a rate of — 0.41"yr^^. Three years later Bosh obtained 
a precession rate of — 0.52"yr~^ when combining the pole position known at that time (1994) with ring plane 



crossings (Bosh et al. 1997). Nicholson & French (1997) and Nicholson et al. (1999) have computed the 
precession rate from 22 reported times of ring plane crossings, reporting a value of — 0.51"yr~^. These values 
refiect the precession rate at the time of the observations, not the long term average. The slow variations in 
Titan's inclination change the torque exerted on Saturn, with a period of about 700 years (Nicholson et al., 
1999). As a result, the measurements must be extrapolated with a model to estimate the long term average 
of Saturn's precession. Currently, the torque seems to be at its minimum value, resulting in a minimum in 



the rate of Saturn's pole precession, ~ 68% of the long-term (secular) value (Vienne & Duriez 1992 1992; 
Nicholson et al., 1999). 

In this section we apply the moment of inertia derived from the interior models, for each assumed rotation 
rate, to derive the long-term precession rate of Saturn's pole. The predicted precession period of Saturn's 
pole due to the solar torque acting on the angular momentum of the Saturnian system is ~ 1.76 x 10^ 



years (French et al. 1993). Recently, Jacobson] (2007) has computed the precession rate from the rigid body 
rotational equations of motion, including the torques from the Sun, Titan and lapetus. Jacobson obtained 
a long term (average value) of — 0.732"yr~^, suggesting a precession period of ^ 1.77 x 10^ years. 



Following the equations in French et al. (1993) we define the precession period by (Ward 19751: 



P 



3J2n'^cose 



(11) 
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where 7 is the planet's moment of inertia, w is Saturn's angular velocity, ns is Saturn's mean motion and e is 
its obliquity. To include the influence of the equatorial satellites, an effective second gravitational moment, 
J2, and an effective moment of inertia, 7', of the Saturian system are used: 



where Mg is Saturn's mass, Rg its equatorial radius, and rrij, aj are the mass and semimajor axis of the j'th 
satellite, respectively. The effective moment of inertia is given by. 



where rij is the satellite's mean motion. The physical parameters used in our computation are summarized 
in Table 5. Saturn's equatorial satellites' masses and radii are given in Table 6. 

[Table. 5] 
[Table. 6] 

To find the precession rate we use the normalized moment of inertia 7 obtained from the empirical interior 
models of section 2. Table 7 presents the computed precession rate of Saturn for different values of 7 and 
rotation period. As previously, the rotation periods range from 10h:32m:35s to 10h:41m:35s. Figure 4. shows 
the normalized moment of inertia and the calculated precession rate values as a function of rotation period. 

[Table. 7] 
[Fig. 4] 

The calculated values can be compared to different long-term precession rates available in literature; 
— 0.7427"yr~^ (French et al., 1993, after modifying the obliquity to the value presented in Table 1), — 0.75"yr~-'^ 
(W'fird & Hamilton, 2004) and — 0.732"yr^^ (Jacobson, 2007). We find that a precession rate of about 
— 0.754"yr^^ is predicted for self-consistent values of rotation rate and derived moment of inertia for all the 
rotation periods considered here. Our models suggest that the long-term value of Saturn's pole precession 
period is ~ 1.72 x 10^ years. 

4.1 Discussion and Conclusions 

We present new models of Saturn's interior with rotation periods between 10h:32m:35s and 10h:41m:35s. 
The models are derived using the 'theory of figures' (Zharkov & Trubitsyn, 1978) with density profiles that 
are represented by a 6th degree polynomial. The interior models fit both the measured Saturnian gravi- 
tational field and the atmospheric model of Fodders & Fegley (1998), providing an empirical polynomial 
density distribution of Saturn's interior (empirical EOS). 




1 



(12) 




(13) 
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Using an EOS of a hydrogen and helium mixture, we find the hydrogen-helium mixing ratio that can best 
match the empirical p(p) relation. Due to uncertainties in the best fit composition value and the empirical 
EOS we present a range of the 'best fit composition'. Since only hydrogen and helium are considered, high-Z 
material is effectively included in the helium mass fraction. When the helium mass fraction exceeds the 
proto-Sun value we evaluate the mass of heavy elements in the interior. The global hydrogen mass fraction 
is found for surface temperatures that range from 130 to 145 K, and rotation periods between 10h:32m:35s 
and 10h:41m:35s. We find that the maximum mass of heavy elements ranges from ~ 6 to 20 M^. 
We look for the 'best fit composition' in different pressure regions. Helium is found to be depleted in the 



upper envelope, in agreement with observations and theoretical models (Conrath & Gautier 2000 Saumon 



& Guillot 


2004 


Guillot 


2005 



interior, the helium mass fraction was found to increase considerably, suggesting that this region is signifi- 
cantly enriched with heavier elements, possibly in the form of a heavy element core. The depletion of helium 
in the upper atmosphere supports the idea of helium differentiation from molecular hydrogen in Saturn's 
envelope (Stevenson, 1975). This process not only explains measurements but it also provides an additional 



energy source insuring that Saturn's evolution is consistent with the age of the solar system (Fortney & 



Hubbard 20031 



Several simplifications have been made in this work. First, the empirical density distribution is given by 
a 6th degree polynomial in fractional radius. As a result, discontinuities in density are smoothed out (see 
Anderson & Schubert, 2007 for further details) and the density in the deep interior might be underestimated, 
especially in the core region. Interior models that include discontinuities (a core) would require a smaller 
enhancement of helium (and any other high-Z material) in the deep interior. However, helium would still be 



expected to be depleted in the outer region and present in larger concentrations in the deep interior (Fortney 



& Hubbard 20031. In addition, it would be desirable to include the high-Z material when calculating the 
'best fit composition' so the distribution of heavy elements within the planet could be better estimated. The 
presence of a heavy element core can then be included as well. 

The estimated hydrogen to helium mixing ratio was found using an adiabatic hydrogen and helium EOS 
based on the calculations of Saumon et al. (1995). Describing the interior by an adiabat is valid as long 
as the planet is fully convective. Convection results in a small superadiabaticity, so the specific entropy is 



expected to be constant throughout the entire planet (Hubbard 1973 Saumon & Guillot 20041. However, 



convection might be suppressed in certain regions due, for example, to the magnetic field, rotation, and 



compositional gradients ( Guillot et al. , 2004 Saumon & Guillot 2004 1 . In the equation of state used in 



this work, the transition from molecular to metallic hydrogen is assumed to be continuous (Saumon et al., 
1995). Thus, if the transition is first order, a discontinuity in entropy can occur, and the assumption of a 
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fully iscntropic planet breaks down (Fortney 20071. Convection can also become inefficient as a result of 
helium differentiation. In this case the entropy of the upper atmosphere will no longer represent the entropy 



of the deep interior, and a two-layer model would be required (Fortney & Hubbard 20031. Finally, different 
equations of state predict densities that can vary by up to 20% at temperatures and pressures relevant to 



the planetary interior ( [Saumon fc Guillot[ |2004[ [Guillot et aL| |2004[ [Militzer et aH|2006 l. 

Improved modeling would be possible when updated atmospheric data for Saturn become available. Updated 
data from Cassini could lead to a smaller error in the atmospheric density resulting in a more accurate 
atmospheric model. In a similar way, a better determination of the temperature at the 1 bar level would 
provide a stronger constraint for the possible adiabats, and therefore, interior models. 

We compute the (long-term) precession rate of Saturn's pole based on the moment of inertia values from the 
empirical interior models. The computed normalized moment of inertia varies with the assumed rotation 
period. We consider rotation periods ranging from 10h:32m:35s to 10h:41m:35s. For the four assumed 
rotation periods we get a precession rate value of ^ — 0.754" yr^^, decreasing slowly with rotation period, 
suggesting that the precession period of Saturn's pole is about ~ 1.72 x 10® years. 
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Saturn Rotation Period 10h:32ni:35s 10h:35m:35s 10h:38m:35s 10h:41m:35s 



a (km) 


60356.2 


60304.3 


60253.4 


60203.5 


R (km) 


58255.4 


58223.4 


58191.9 


58161.1 


q 


0.158851 


0.156949 


0.155085 


0.153257 


m 


0.142835 


0.141256 


0.139706 


0.138182 


Po (kg m"3) 


686.276 


687.409 


688.525 


689.621 


Po (Mbar) 


4.46847 


4.47832 


4.48800 


4.49754 


J2 (io-«) 


16276.6 


16304.6 


16332.2 


16359.2 


Ji (10-6) 


-934.2 


-937.4 


-940.6 


-943.7 


Je (10-') 


85.9 


86.4 


86.8 


87.2 



Table 1: Fifth order (in smallness parameter, q defined in text) reference geoid for four rotation periods 
that span the six-minute interval of possible periods. The polar radius is fixed at the value 54,438 km, the 



polar radius of the 100 mbar isosurface (Lindal et al. 19851. The harmonic coefficients J2„ are obtained 



from the measured values which are given for a reference equatorial radius of 60,330 km (Jacobson et 
al. 2006), according to a^" J2„=(60,330 km)^" J2„(measured). Since the equatorial radius, a of the reference 
geoid depends on the rotation period, so do the values of J2n- The equatorial radius, a and the values of 
J2n are determined iteratively using the above relation and a second equation derived from equation (1) 
evaluated at the pole and the equator. The values of J271 from Jacobson et al. (2006) in units of 10^' are 
J2 = 16290.71 ± 0.27, J4 = —935.8 ± 2.8, Jq = 86.1 ± 9.6. The parameter po is a characteristic pressure 
defined in the text. 
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Saturn Rotation Period 10h:32ni:35s 10h:35m:35s 10h:38m:35s 10h:41m:35s 



ko 






6.459802 


6.583734 


6.712154 


6.845172 


k2 






-35.442449 


-38.344462 


-41.325925 


-44.380760 


ka 






12.828048 


18.625924 


24.583873 


30.657281 


k4 






110.378345 


110.132162 


109.835175 


109.590534 


k5 






-156.131034 


-162.169863 


-168.296182 


-174.607977 


AJ2 


(10- 


-6) 


0.05 


0.06 


0.06 


0.07 


AJ4 


(10- 


-6) 


0.88 


0.95 


1.02 


1.09 


AJg 


(10- 


-6) 


3.47 


3.75 


4.03 


4.32 


7 






0.222711 


0.223809 


0.224902 


0.225990 



Table 2: Best-fit values for the polynomial coefficients that describe the variation of normalized density 



77 as a function of normalized mean radius /3. The coefficients fit a model atmosphere (Lodders & Fegley 



1998) between 1000 mbar and 100 mbar with an assumed error of 100% on the atmospheric density. They 
also fit the gravitational harmonics J2, J a, and Jg with an assumed error given by the covariance matrix for 
Jin (Jacobson et al. 2006). Saturn's normalized axial moment of inertia 7 ~ C jMR^ is computed from the 
density polynomial. The converged residuals for each rotation period are given by AJ2„. The fit to both 
the atmosphere and the harmonics is best for the shortest period of 10h:32m:35s, but all fits are acceptable, 
within one standard error of the measured value. The residuals in the normalized atmospheric density are 
on the order of 10^^ and are not shown. All four models fit the atmospheric density to well within the 
assumed error of 100%. The failure to fit the harmonics exactly is a measure of the incompatibility between 
the model atmosphere and the gravitational field. 



Saturn Rotation Period 


T(l bar) [K] 


Best Fit Value 


Composition Range 


Maximum Mass of Heavy Elements [Mgj] 




130 


0.82 


0.66-0.98 


6.2 




135 


0.74 


0.55-0.93 


16.6 


10il:32m:35s 


140 


0.71 




18.6 


0.53-0.89 




145 


0.67 


0.50-0.84 


21.4 




130 


0.80 


0.65-0.95 


7.1 




135 


0.74 


0.57-0.91 


14.8 


10h:35m:35s 


140 


0.73 




12.0 


0.60-0.86 




145 


0.67 


0.52-0.82 


19.5 




130 


0.76 


0.59-0.93 


12.9 




135 


0.74 


0.59-0.89 


12.9 


10h:38m:35s 


140 


0.71 




13.8 


0.58-0.84 




145 


0.66 


0.52-0.80 


19.5 




130 


0.75 


0.58-0.92 


13.8 




135 


0.73 


0.59-0.87 


12.9 


10h:41m:35s 


140 


0.67 




19.5 


0.52-0.82 




145 


0.65 


0.51-0.79 


20.5 



Table 3: The best fit of the hydrogen mass fraction for diflterent rotation rates of the planet's interior and temperatures 
at the 1 bar level ranging from 130 to 145 K. 
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Temperature (1 bar) 


Rotation Period 


best fit X value 


best fit X value 


best fit X value 






logP{Mbar) = -3 ~ 1.12 


logP{Mbar) =(-6)-(-4) 


logP{Mbar) = ~ 1.12 


130 K 


10h:32m:35s 


0.58 ± 0.16 


1.00 ± 0.02 


0.28 ± 0.05 


135 K 


10h:32m:35s 


0.54 ± 0.17 


0.99 ± 0.03 


0.27 ± 0.06 


140 K 


10h:32m:35s 


0.53 ± 0.16 


0.97 ± 0.05 


0.27 ± 0.06 


145 K 


10h:32m:35s 


0.53 ± 0.14 


0.93 ± 0.05 


0.27 ± 0.05 


130 K 


10h:35m:35s 


0.57 ± 0.14 


1.00 ± 0.02 


0.28 ± 0.06 


135 K 


10h:35ni:35s 


0.54 ± 0.13 


1.00 ± 0.02 


0.27 ± 0.05 


140 K 


10h:35m:35s 


0.53 ± 0.13 


0.98 ± 0.05 


0.27 ± 0.05 


145 K 


10h:35m:35s 


0.51 ± 0.13 


0.91 ± 0.05 


0.27 ± 0.05 


130 K 


10h:38m:35s 


0.54 ± 0.14 


1.00 ± 0.02 


0.28 ± 0.04 


135 K 


10h:38m:35s 


0.53 ± 0.14 


1.00 ± 0.02 


0.27 ± 0.05 


140 K 


10h:38m:35s 


0.53 ± 0.13 


0.96 ± 0.04 


0.27 ± 0.05 


145 K 


10h:38m:35s 


0.51 ± 0.15 


0.89 ± 0.04 


0.27 ± 0.05 


130 K 


10h:41m:35s 


0.54 ± 0.14 


1.00 ± 0.03 


0.28 ± 0.07 


135 K 


10h:41ni:35s 


0.53 ± 0.13 


1.00 ± 0.03 


0.28 ± 0.07 


140 K 


10h:41m:35s 


0.51 ± 0.13 


0.94 ± 0.05 


0.27 ± 0.06 


145 K 


10h:41m:35s 


0.48 ± 0.13 


0.87 ± 0.04 


0.27 ± 0.06 



Table 4: Best fit of composition for different pressure regions assuming different surface temperatures and rotation 
periods. 



GMsaturn (Ws-^) 


37,931,207.7(1) 


GM5™ (km3s-2) 


1.3712440018x1020(2) 


h 


16290.71(3) 


Rs, Saturn equatorial radius (km) 


60268(2) 


as, Saturn semimajor axis (km) 


1.460268x10^(2) 


e, obliquity (deg) 


26.73919('i) 



Table 5: Pliysical parameters: (^'Anderson & Scfiubert (2007), '^'JPL data,'^' Jacobson et al. (2006),'"*' Jacobson 
(2007) 



Satellite 


GM km3s-2 


Semimajor axis (lO^km) 


Mimas 


2.530 


185.54 


Enceladus 


7.210 


238.04 


Tethys 


41.210 


294.67 


Dione 


73.113 


377.42 


Rhea 


154.07 


527.07 


Titan 


8978.19 


1221.87 


Hypreion 


0.37 


1500.88 


lapetus 


120.50 


3560.84 



Table 6: Satellite data, JPL database: http://ssd.jpl.nasa.gov 
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Rotation Ptaiod 




(\)iiil)ut(Hl Pr(K:(»hioii Rate 


10h:32m:35s 


0.222711 


-0.7544"yr-i 


10h:35m:35s 


0.223809 


-0.7543"yr~^ 


10h:38m:35s 


0.224902 


-0.7541"yr-i 


10h:41m:35s 


0.225990 


-0.7540"yr-i 



Table 7; Precession rates for different values of normalized moment of inertia and rotation periods. 



Error of the p(p) relation (Rotation Period- 10 h, 32 m, 35 s) 




-6 -5 -4 -3 -2 -1 

Log p (Mbar) 



Figure 1: Saturn's empirical EOS with the error included. The solid line shows the computed p{p), and the 
dotted and dashed-dotted curves present the p{p) relation when the systematic error is added and subtracted, 
respectively. The area between these two curves represents the empirical EOS from the interior model. The 
interior model presented refers to a rotation period of 10h;32m:35s. 
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-2.0 -1.5 -1.0 -0.5 0.0 0.0 0.2 0.4 0.6 0.8 1.0 

Log p (Mbar) Log p (Mbar) 



Figure 2: Saturn's empirical p{p) relation. The solid and dotted curves represent rotation periods of 
10h:41m:35s and 101i:32m:35s, respectively. 
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Rotation Period- 10 h, 32 m, 35 s Rotation Period- 10 ti, 35 m, 35 s 




-6 -5 -4 -3 -2-10 1 -6 -5 -4 -3 -2-10 1 

Log p (Mbar) Log p (Mbar) 



Rotation Period- 10 h, 38 m, 35 s Rotation Period- 10 h, 41 m, 35 s 




-6 -5 -4 -3 -2-10 1 -6 -5 -4 -3 -2-10 1 

Log p (Mbar) Log p (Mbar) 



Figure 3: Saturn's empirical EOS for rotation periods ranging from 10h:32m:35s to 10h:42m:35s for dif- 
ferent surface temperatures. Each plot shows the difference between the empirical EOS and the physical 
one, ALogp = {Logpn/He — LogpempiHcai)- The dotted, solid, dashed-dot and dashed curves represent 
temperatures of 130, 135, 140 and 145 K at the 1 bar level, respectively. 
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10.54 10.56 10.58 10.60 10.62 10.64 10.66 10.68 

Rotation Period (hr) 

Figure 4: Saturn's normalized moment of inertia 7 (dashed) and precession rate (solid) as a function of 
rotation period. 
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